##-----------------------------------------##
## Lauren Peritz
## ISQ 2019 - WTO Compliance
## Replication File Part 4
## Use "Peritz_ISQ2019_WTO4.dta"              
## Update: 2019.11.19                                 
###########################################

#################################################################################
#### 1: Preliminaries -----------------------------------------------------------

rm(list=ls(all=TRUE))

# Load packages
library(Synth)
library(reshape)
library(foreign)
library(scales)
library(xtable)
library(haven)


# Set working directory
setwd("~/Dropbox/WTO_Compliance/Replication")

# Import data
dat <- as.data.frame(read_dta("~/Dropbox/WTO_Compliance/Replication/Peritz_ISQ2019_WTO4.dta"))


#################################################################################
#### 2: Prepare SCM Calculations ------------------------------------------------


# Identify disputes
ds.list <- levels(as.factor(dat$Dispute))

# Custom color for plots
tr.gray <- alpha("gray", alpha=.3)

# Default predictors (covariates) to use in SCM
pred.default <- c("gdp","gdp_growth","gdppc","agriculture_value_added",  
                  "industry_value_added","manufacturing_value_added",
                  "services_value_added","trade","unemployment_t")



#################################################################################
#### 3: Functions ------------------------------------------------


## scm.analysis()

scm.analysis <- function (DS.No, Resp, dat, predictors.list, window = 5){
  
  # Parameters
  DS <- dat[dat$Dispute==DS.No,]
  
  year.min <- min(DS$Year)
  year.max <- min((mean(DS$RPTExpireYear)+window), max(DS$Year))
  year.rpt <- mean(DS$RPTExpireYear)
  donor.pool <- unique(DS$PartnerCode[DS$PartnerCode!=Resp])
  
  # subset Dataset
  DS <- DS[DS$Year<=year.max,]
  attach(DS)
  
  # Compute Log-Odds as log(Share.Respondent/Share.Next.Largest)
  DS.sub <- DS[DS$PartnerCode!=Resp,c(2,5,11)]
  control.means <- by(DS.sub$Share, INDICES=as.factor(DS.sub$PartnerCode), FUN = mean)
  denom <<- DS.sub[DS.sub$PartnerCode==labels(control.means[control.means==max(control.means)]),]
  denom.long <- rep(denom$Share,length(levels(as.factor(DS$PartnerCode))))
  DS$l.odds.rat <- log(DS$Share/denom.long)
  
  # SCM calculation
  dataprep.out <<- dataprep(foo=DS, 
                            predictors=predictors.list,
                            dependent = "l.odds.rat", 
                            unit.variable = "PartnerCode", time.variable = "Year",
                            treatment.identifier = Resp, 
                            controls.identifier = donor.pool,
                            time.predictors.prior = c(year.min:year.rpt),   
                            time.optimize.ssr = c(year.min:year.rpt),
                            time.plot = (year.min:year.max),
                            unit.names.variable = "PartnerName")
  synth.out <<- synth(dataprep.out)

  # Summary tables 
  synth.tables <<- synth.tab(dataprep.res = dataprep.out,synth.res = synth.out)
  print(synth.tables)

}

# ## ds.stat()
# ## Calculates compliance statistic, transformed to trade share

ds.stat<-function(DS.No){
  	all <-length(dataprep.out$tag$time.plot)
  	c <-length(dataprep.out$tag$time.predictors.prior)
	  gaps <- exp(dataprep.out$Y1plot)*denom$Share-exp(dataprep.out$Y0plot%*%synth.out$solution.w)*denom$Share
  	gaps.diff<-round(mean(gaps[c:all,])-mean(gaps[1:c,]),6)
  	gaps.sd <- (round( sd(as.vector(gaps[1:c,]))  ,6))/2
  	stat <- as.data.frame(cbind(DS.No, gaps.diff, gaps.sd))
  	return(stat)
	}# end ds.stat()


	
## path.plot.exp() 
## Plots SCM, transformed back to trade share

path.plot.exp<-function (synth.res = NA, dataprep.res = NA, tr.intake = NA, 
    Ylab = c("Y Axis"), Xlab = c("Time"), Ylim = NA, Legend = c("Treated", 
        "Synthetic"), Legend.position = c("topright"), Main = NA, 
    Z.plot = FALSE, cex.l=6/7) 
{
    if (Z.plot == FALSE) {
        if (sum(is.na(dataprep.res$Y1plot)) > 0) {
            stop("\n\n#####################################################", 
                "\nYou have missing Y data for the treated!\n\n")
        }
        if (sum(is.na(dataprep.res$Y0plot)) > 0) {
            stop("\n\n#####################################################", 
                "\nYou have missing Y data for the controls!\n\n")
        }
        y0plot1 <- dataprep.res$Y0plot %*% synth.res$solution.w
        if (sum(is.na(Ylim)) > 0) {
            Y.max <- max(c(y0plot1, dataprep.res$Y1plot))
            Y.min <- min(c(y0plot1, dataprep.res$Y1plot))
            Ylim <- c((Y.min - 0.3 * Y.min), (0.3 * Y.max + Y.max))
        }
        plot(dataprep.res$tag$time.plot, exp(dataprep.res$Y1plot)*denom$Share, 
            t = "l", col = "black", lwd = 2, main = Main, ylab = Ylab, 
            xlab = Xlab, xaxs = "i", yaxs = "i", ylim = Ylim)
        lines(dataprep.res$tag$time.plot, exp(y0plot1)*denom$Share, col = "black", 
            lty = "dashed", lwd = 2, cex = 4/5)
    }
    else {
        z0plot <- dataprep.res$Z0 %*% synth.res$solution.w
        if (sum(is.na(Ylim)) > 0) {
            Y.max <- max(c(z0plot, dataprep.res$Z1))
            Y.min <- min(c(z0plot, dataprep.res$Z1))
            Ylim <- c((Y.min - 0.3 * Y.min), (0.3 * Y.max + Y.max))
        }
        plot(dataprep.res$tag$time.optimize.ssr, exp(z0plot)*denom$Share, t = "l", 
            col = "black", lwd = 2, main = Main, ylab = Ylab, 
            xlab = Xlab, xaxs = "i", yaxs = "i", ylim = Ylim)
        lines(dataprep.res$tag$time.optimize.ssr, exp(dataprep.res$Z1)*denom$Share, 
            col = "black", lty = "dashed", lwd = 2, cex = 4/5)
    }
    abline(v = tr.intake, lty = 3, col = "black", lwd = 2)
    if (sum(is.na(Legend)) == 0) {
        legend(Legend.position, legend = Legend, lty = 1:2, col = c("black", 
            "black"), lwd = c(2, 2), cex = cex.l, bty="n")
    }
}



#################################################################################
#### 4: SCM for Disputes 8 and 75 (FIGURE 2) --------------------



## DS8: EU v. Japan - Alcoholic Beverages	
scm.analysis(DS.No=8, Resp=392, dat=dat, predictors.list=c(pred.default,"polity"))
ds.stat(DS.No=8)[2:3]
detach(DS)

## Plot FIG2b_ds8
par(mar=c(5,5,6,2))
path.plot.exp(dataprep.res = dataprep.out,synth.res = synth.out, 
       Ylim = c(-.015, .10),
	   Ylab = "EU Alcohol Exports to Japan (Share)",
       Xlab = "Year",
       Legend = c("Japan", "Synthetic Japan"))
rect(xleft=1995+5/12, xright=1998+1/12,ybottom=-0.015, ytop=.10, col=tr.gray, border = tr.gray)
abline(v=1998+1/12, lwd=2, col="grey35")
mtext("EU v. Japan - Alcoholic Beverages (DS8)", side=3, line=1, cex=1.4)



## DS75: EU v. Korea - Alcoholic Beverages  
scm.analysis(DS.No=75, Resp=410, dat=dat, predictors.list=c(pred.default,"polity"))
ds.stat(DS.No=75)[2:3]
detach(DS)

## Plot FIG2a_ds75
 par(mar=c(5,5,6,2))
 path.plot.exp(dataprep.res = dataprep.out,synth.res = synth.out, 
               Ylim = c(0, .12),
               Ylab = "EU Alcohol Exports to Korea (Share)",
               Xlab = "Year",
               Legend = c("Korea", "Synthetic Korea"),
               Legend.position = "topleft")
 rect(xleft=1997+3/12, xright=2000+1/12,ybottom=0, ytop=.12, col=tr.gray, border = tr.gray)
 abline(v=2000+1/12, lwd=2, col="grey35")
 mtext("EU v. Korea - Alcoholic Beverages (DS75)", side=3, line=1, cex=1.2)


# End